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Abstract 


The square-lattice quantum Heisenberg antiferromagnet displays a pronounced anomaly 
of unknown origin in its magnetic excitation spectrum. The anomaly manifests itself only 
for short wavelength excitations propagating along the direction connecting nearest neigh¬ 
bors. Using polarized neutron spectroscopy, we have fully characterized the magnetic fluc¬ 
tuations in the model metal-organic compound CFTD, revealing an isotropic continuum 
at the anomaly indicative of fractional excitations. A theoretical framework based on the 
Gutzwiller projection method is developed to explain the origin of the continuum at the 
anomaly. This indicates that the anomaly arises from deconfined fractional spin-1/2 quasi¬ 
particle pairs, the 2D analog of ID spinons. Away from the anomaly the conventional 
spin-wave spectrum is recovered as pairs of fractional quasiparticles bind to form spin-1 
magnons. Our results therefore establish the existence of fractional quasiparticles in the 
simplest model two dimensional antiferromagnet even in the absence of frustration. 


A fascinating manifestation of quantum mechanics in the presence of many-body interac¬ 
tions is the emergence of elementary excitations carrying fractional quantum numbers. Frac¬ 
tional excitations were a central ingredient to understand the anomalous quantum Hall effect [Q]|, 
and have been investigated in a range of systems including conducting polymers A3, bilayer 
graphene A3], cold atomic gases 0, and low-dimensional quantum magnets A3 ID- Among 
the latter class of systems, the spin- 1/2 chain quantum Heisenberg antiferromagnet (HAF) is 
perhaps the simplest realization of this paradigm, for which the ground-state and the excitation 
spectrum known exactly AZl 13 0- The elementary excitations of an antiferromagnet in the large- 
spin, semi-classical case are well-defined spin-waves (magnons) created by AS = 1 processes, 
while in the quantum spin- 1/2 chain a AS = 1 process generate pairs of unbound fractional 
quasi-particles known as spinons each carrying a S = 1/2 quantum number AlEJ. While the ex¬ 
istence of spinons in the spin- 1/2 HAF chain has been established in a number of experimental 
studies on model ID materials [flOlITTl 12], establishing their existence in higher dimensions is 
an ongoing challenge AH. To date the main candidate systems comprise antiferromagnets with 
frustrated geometries, including triangular AT3l or kagome Ifl4l '15, 16:] lattices. In this work, 
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we take a frustration-free route and focus on the spin-1/2 square-lattice quantum Heisenberg 
antiferromagnet (SLQHAF), one of the most fundamental models in magnetism. It is defined 
by the Hamiltonian 



( 1 ) 


(*d> 


where the antiferromagnetic exchange J is restricted to nearest-neighbor spins We pro¬ 
vide experimental and theoretical evidence that even in this simplest of 2D models deconfined 
fractional S — 1/2 quasiparticles exist and are responsible for the very significant deviations of 
the high-energy, short-wavelength spin dynamics from conventional spin-wave theory. 

It may seem surprising that the SLQHAF is a candidate for hosting fractional excitations 
as at a superficial level its magnetic order and excitations resemble those of a classical system. 
It is well established that at zero temperature quantum fluctuations are insufficient to prevent 
the development of collinear antiferromagnetic order ff71 . The elementary excitations of this 
Neel state, when calculated using semi-classical spin wave theory (SWT), can be described as 
S — 1 bosonic quasiparticles, known as magnons: the one-magnon dispersion is gapless [18, 
ED, with two-magnon excitations occupying a continuum at higher energy. The interaction 
between magnons is relatively weak and leads to an upward renormalization of the magnon 
energy [20, 21 ] and to scattering between two-magnon states Il22t[23 |. 

While none of the above properties suggest the existence of quasiparticle fractionalization, 
quantum effects are nevertheless far from negligible in the SLQHAF. This is evidenced by the 
observation that quantum zero-point fluctuations reduce the staggered moment to only 62% of 
its fully ordered value S — 1/2 [24, 25]. This suggests that the SLQHAF may in fact be close 
to a state lacking spin-symmetry breaking, such as the resonating-valence-bond (RVB) state 
proposed for the cuprate realization of this model by Anderson fl26ll . In particular, fractional 
spin excitations present in the RVB state may be relevant for the spin dynamics in the Neel state, 
especially at high energies. Indeed, analytical theories using bosonic Il27ll or fermionic If28ll29 1 
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Figure 1: (a) Overview of the magnetic excitation spectrum of CFTD. Momentum and en¬ 

ergy dependence of the (total) dynamic structure factor <S(q, u) measured by time-of-flight 
inelastic neutron scattering. Square boxes (black dashed) highlight the (jr, 0) and (7t/2,7t/2) 
wave-vectors, (b) and (d) Corresponding distributions of real-space fractional quasiparticle- 
pair separations, as calculated in the |SF) variational state (Eq. [6]), evidencing respectively the 
unbound and bound nature of the pair excitations, (c) Pictorial representation of a quasiparticle 
pair excitation, (e) Pictorial representation of a spin-wave excitation (magnon). 


fractional quasiparticles have long been proposed. By analogy with the ID case, these are 
referred to as spinons. 

The magnetic excitation spectrum of various realizations of the SLQHAF have been investi¬ 
gated using neutron spectroscopy, including parent compounds of the high-T c cuprates such as 
La 2 Cu0 4 L3QU3T] and the metal-organic compound Cu(DC00) 2 4D 2 0 (CFTD) [32, 33] con¬ 
sidered here. These experiments have established that while SWT gives an excellent account 
of the low-energy spectrum, a glaring anomaly is present at high energy for wavevectors in 
the vicinity of (n, 0), where q = ( q x ,Q y ) is expressed in the square-lattice Brillouin-zone of 
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unit-length 2 tt. The anomaly is evident as a dramatic wipeout of intensity of the otherwise 
sharp excitations, Fig. [7}i for CFTD and Ref. PTll for La 2 Cu0 4 , and as a 7% downward disper¬ 
sion along the magnetic zone-boundary connecting the (7 t/2,7t/ 2) and (zr, 0) wavevectors for 
CFTD (see also Fig. Iff)- Unambiguously identifying the origin of this effect is complicated 
by the presence, in some materials, of small additional exchange terms such as electronic ring- 
exchange or further neighbor exchange. In contrast, the deviations observed in CFTD agree 
with numerical results obtained by series expansion 0411351 . quantum Monte Carlo [36. 37J, 
and exact diagonalization If38ll methods for the model of Eq. (|T|), proving that the anomaly is in 
this case intrinsic ll33l . Given the fundamental nature of the SLQHAF, it is clearly desirable to 
develop a microscopic understanding of the origin of the anomaly. 

Here we present polarized neutron scattering results on CFTD which establish the existence 
of a spin-isotropic continuum at (7r, 0), which contrasts sharply with the dominantly longitu¬ 
dinal continuum at (tt/ 2, tt / 2) and with expectations for a magnetically ordered ground state. 
Using a fermionic description of the spin dynamics based on a Gutzwiller-projected variational 
approach, we argue that this is a signature of unbound pairs of fractional S — 1/2 quasiparticles 
(Fig.|TJ> and[TJ:) at the (tt, 0) wave vector. At other wave vectors, including (tt/2, it/2) (Fig. Iff), 
our approach yields bound pairs these fractional quasiparticles and so recovers a conventional 
magnon spectrum, in agreement with SWT (Fig.|T] 2 ). 

Neutron scattering experiments were performed on single crystals of CFTD using unpolar¬ 
ized time-of-flight spectroscopy (Fig. |Tj) and triple-axis spectroscopy with longitudinal polar¬ 
ization analysis (see Supplementary Materials). The results of our polarized experiment are 
presented in Fig. [2] through the energy dependence of the diagonal components of the dynamic 
structure factor. By combining wave-vectors from different equivalent Brillouin zones (see 
Supplementary Materials), we can reconstruct the total dynamic structure factor (Fig. [2 ^ and 
e) and separate contributions from spin fluctuations that are transverse to and along (Fig. @>-c 
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Figure 2: Summary of the polarized neutron scattering data. Energy dependence of the to¬ 
tal, transverse and longitudinal contributions to the dynamic structure factor at constant wave 
vectors (a-d) q = ( 7 r, 0) and (e-h) q = (7t/2, 7t/ 2) measured by polarized neutron scattering 
on CFTD. The solid lines indicate resolution-limited Gaussian fits, while the dashed lines are 
empirical lineshapes used as guides-to-the-eye. 


and f-e) the ordered moment. Within SWT, the resulting transverse and longitudinal spectra are 
dominated by one-magnon and two-magnon excitations, respectively. At (7 t/2 , 7 t/2 ) and at an 
energy of u/ J = 2.38(2) we observe a sharp, energy resolution limited peak (Acu = 1.47(5) meV 


6 


































= 0.24(1) J, FWHM) which is the signature of a long-lived, single-particle excitation (Fig.[2]e). 
Most of the observed spectral weight is in the resolution-limited peak of the transverse channel 
5 ± (q, u) = S xx (q,uj) + S yy ( q,cu) (Fig.[2]f), while a weak continuum extends from u/ Jm 2.3 
to 3.4 with a maximum around ce/ J ~ 2.6 in the longitudinal channel, <S 2Z (q, oS) (Fig. [2^). In 
contrast, the response at (n, 0) displays a pronounced high-energy tail, starting right above the 
peak maximum at c of J = 2.19(2) and extending up to to/J «3.8. This tail carries 40(12)% of 
the total spectral weight at ( 7 r, 0) (Fig. [2J1) and is evident in both the transverse (Fig. |2]b) and 
longitudinal (Fig. |2]c) channels. To isolate the continuous component in the transverse chan¬ 
nel we subtract resolution-limited Gaussians corresponding to sharp, single-particle responses, 
with the results shown in Fig. [2]cl and h. This analysis reveals the important fact that the trans¬ 
verse continuum at (7 r, 0) is within error twice the longitudinal contribution (Fig.[2)i). Thus we 
can conclude that the continuum at ( 7 r, 0) arises from correlations which are isotropic in spin 
space with S-^q, u) = 2S zz (q, ui), while by contrast the continuum contribution at (7r/2, 7t/2) 
is fully contained in the longitudinal channel (Fig.[2ji). 

The pronounced asymmetric and non-Lorentzian lineshape of the continuum at ( 7 r, 0) can¬ 
not be accounted for by conventional effects, even including instrumental resolution. SWT 
predicts that magnon interaction transfer up to 20% of the transverse spectral weight at the 
zone-boundary from the sharp one-magnon peak to a higher energy continuum of three-magnon 
states Ii23ft . However, the resulting lineshape differs radically from our observations, does not 
coincide with the longitudinal response, and does not appear to depend significantly on the 
wave vector along the zone boundary. Spontaneous magnon decays can in principle produce 
an asymmetric lineshape but are prohibited in this case by the collinearity of the magnetic or¬ 
der [|22l [391 - Instead, recent quantum Monte Carlo work [|40l suggests to look for explanations 
of the continuum contribution to the dynamic structure factor at ( 7 r, 0) involving the deconfine¬ 
ment of fractional excitations. This is further motivated by the observed coexistence of sharp 
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two-spinon bound states with a broad multi-spinon continuum, at comparable energy ranges but 
different wave-vectors, in the quasi-2D materials CS 2 CUCI 4 [13, 41] and LiCuVOj |f42] made 
of strongly-coupled Heisenberg chains. 

In order to explore whether fractionalization of magnons can account for the ( 7 r, 0) anomaly 
in the SLQHAF, we use a theoretical approach based on Gutzwiller-projected variational wave 
functions[43, 44]. In this approach, spin operators are transformed into pairs of S = 1/2 
fermionic operators so that Eq.|T]becomes 



( 2 ) 


constant, 


where c\ a () creates (annihilates) an electron with spin a at site 1 . This transformation em¬ 


beds the original spin Hilbert space into an electronic Hilbert space which also contains non¬ 
magnetic sites occupied by zero or two electrons. As a result, Eqs[j]and[2]are only equivalent on 
the restricted electronic subspace with half electron filling and no empty sites or double occu¬ 
pancies. This constraint can be enforced exactly by the so-called Gutzwiller projector P G . The 
advantage of this approach is that pairs of fractional S — 1/2 quasiparticles (for the original spin 
Hamiltonian) can be naturally constructed as particle-hole excitations in the electronic space, 
projected a-posteriori by P G onto spin configurations with exactly one electron per site. 

The quartic electronic operator in Eq.[2]is treated by a mean-field decoupling where the av¬ 
erages (P ia c i(T ) and {c\ a Cj a ) are considered. We adopt the following Ansatze for their real-space 
dependencies : (, cl a c i(T ) corresponds to a staggered Neel order parameter (N) and (c| CT c JCT ) to a 
staggered flux (SF) threading square plaquettes of the lattice [[45, [46; 4711 (see Supplementary 
Materials for exact definitions and more details). To each average corresponds a variational 
parameter whose value is optimized to minimize the energy (Eq. |T]) of the Gutzwiller-projected 
state, |SF+N) = / //T’mf)- The corresponding mean-field electronic ground-state Amf) con- 
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tains empty and doubly occupied sites and reads 


IV’mf) = JJ 7kt-7i4-|0), (3) 

keMBZ 

where | 0 ) is the electron vacuum and where the 7 operators are linear combinations of c^ a 
operators that diagonalize the mean-field electronic Hamiltonian. The product over the wave- 
vector k is restricted to the magnetic Brillouin zone (MBZ), a result of the antiferromagnetic 
unit-cell-doubling. Consequently “±” denotes the band index. In the present case of half elec¬ 
tron filling, the ” band is fully occupied, and there is a finite gap to the empty “+” band 
for non-zero Neel order-parameter. The overall minimization procedure is carried out numer¬ 
ically using Variational Monte Carlo [1441 and leads to a |SF+N) state with variational energy 
-EgF+N = —0.664J and staggered moment 0.75,5 per site [43, 48]. This can be compared to 
more precise Green’s function Monte Carlo studies for Eq.|T]that obtained — 0.669 J and 0.615,5' 
for the ground-state energy and the staggered moment, respectively [l49l 50 1 . 

The optimized |SF+N) state, while giving a good estimate for the ground-state energy, does 
not have the correct long-distance behaviour for the transverse equal-time correlator ( S + (0) S~ (r)), 
predicted by SWT to decay as a power-law OTTl . Instead, as the excitation spectrum of the 
mean-field electronic ground-state is gapped, (S + (0)S~(r)) decays exponentially after projec¬ 
tion. We conjecture that the asymptotic behavior of the spin correlator is important for the 
deconfinement of fractional excitations. To obtain insight into the influence of long-distance 
spin fluctuations, we consider a distinct variational state, |SF), for which the finite staggered- 
flux is retained but the Neel order is reduced to zero. |SF) is a quantum spin-liquid singlet 
of variational energy E$f = —0.638J that displays a power-law decay of its transverse spin 
correlations |5T1[52I. 

We now turn to the construction of transverse (S = 1) spin excitations for the above varia¬ 
tional states, aiming at comparing their respective dynamic structure factor with the results of 
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Figure 3: A schematic representation of local-spin-flip and spatially separated quasiparticle 
pair excitations in the Gutzwiller-projected approach, (a) The mean-field wave function IV’mf) 
is shown as a resonating-valence-bond liquid (for a better visualization, all singlets are shown 
as nearest-neighbor and the Neel order is ignored). Configurations containing doubly occupied 
sites (right hand side) are discarded by the Gutzwiller projection P G . (b) Local spin flips create 
triplets out of resonating singlets. Configurations from I^mf) originally containing doubly oc¬ 
cupied sites are still projected out (right hand side), (c) Nonlocal quasiparticle-pair excitations 
are constructed as projected particle-hole excitations. At a nonzero separation r, they contribute 
by annihilating a doubly occupied site with a hole, leaving two separated spins-j\ After projec¬ 
tion, the only configurations left are the ones constructed from IV’mf) that contained one empty 
and one doubly occupied site (right hand side). 

Fig. [2} The variational transverse spin excitations are obtained as superpositions of projected 
particle-hole excitations with momentum q: 

|q,™,+)= ^kql^q)’ |k,q) = PcTit+Tk-qi-I^MF) (4) 

keMBZ 

where the states |k, q) are generated by destroying a spin-down quasiparticle in the ” band 
and creating a spin-up quasiparticle in the “+” band. The coefficients </>£ q are obtained by di¬ 
agonalizing the original Hamiltonian (Eq. |T|) projected onto the non-orthonormal set of states 
|k, q) and correspond to the eigen-energies E+ n . Expressing the Fourier-space quasiparticle 
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operators 7 ko - ± using the real-space c ia operators, we note that the variational spin excitations 
contain both local spin flips S+P G IV’mf) = / (Fig. [ijb) and spatially separated 

particle-hole excitations, P G (Fig. 3je) . The dynamic structure factor of the trans¬ 

verse spin excitations is calculated as 


S + -(q, w) = £ |{q. n, +|SJ|GS)| 2 <5 (w - £+ + E as ) (5) 

n 

where |GS) stands either for |SF+N) or |SF). We use the identity S = S 1 = S ~ + valid for 
both variational ground-states to compare the transverse dynamic structure factor of the vari¬ 
ational states |SF+N) and |SF) with the experimental results presented in Fig. [2j A similar 
approach also allows to obtain the longitudinal (S = 0) dynamic structure factor (see Supple¬ 
mentary Materials). 

The transverse dynamic structure factor ^(q, u) of the |SF+N) state is shown in Fig.jdja), 
as obtained by variational Monte Carlo on a finite lattice of 24 x 24 sites. The dominant feature 
of the spectrum is a low-energy magnon-like mode, which resembles the experimental results of 
Fig#- In particular, our calculation produces a dispersion along the magnetic zone-boundary 
in better quantitative agreement with the 7% dispersion observed in Ref. 031 than any other the¬ 
oretical method, Figs. [4]: an dp. This confirms that magnons can be quantitatively interpreted 
as bound pairs of fractional S'= 1/2 quasiparticles. 

However the |SF+N) transverse dynamic structure factor exhibits a gap at (n,n) and no 
continuum above the magnon branch at (n, 0). We believe that this is an artifact of replacing 
the spontaneous symmetry breaking by a Neel mean-field order parameter: this ansatz, as men¬ 
tioned above, distorts the long-distance asymptotics of spin correlations. Indeed, if we reduce 
the Neel mean-field parameter of the |SF+N) state, then the high-energy excitations at ( 7 r, 0) 
move down in energy (see Supplementary Materials). At the vanishing Neel field (i.e., in the 
|SF) state) they evolve into a succession of modes distributed on an extended energy range 
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Figure 4: Calculated transverse dynamic structure factor for the |SF+N) (a) and |SF) (b) states 
with lattice sizes of 24 x 24 and 32 x 32 respectively, (c) The magnon-like dispersion extracted 
from a (red points) compared to the experimental CFTD data Il33l (blue squares), spin-wave 
theory with first-order (solid black line) and third-order 1(531 (green triangles) 1/S corrections, 
series expansion lf35l (dashed purple line) and quantum Monte Carlo l(37l (cyan diamonds). 
The experimental data are scaled using J = 6.11 meV. (d) Zoom-in on the magnon-like mode 
dispersion along the magnetic zone-boundary. 


above the spin-wave mode (shown in Fig. B> for a 32 x 32 lattice). This behavior contrasts 
the situation at (zr/2, zr/2) where the high-energy transverse excitations completely lose their 
spectral weight on reducing the Neel field and only the spin-wave mode remains in the |SF) 
state. We therefore suggest that the continuum of excitations observed at (7r, 0) is related to 
power-law transverse spin correlations and that it corresponds to deconfined fractional spin-1/2 
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quasiparticles. 

To support this interpretation, we consider in Fig. [5}r and[5]b the system-size dependence 
of ^(q, oj). While the excitations at (tt/2, tt/2) form a single sharp mode with energy and 
intensity nearly independent of the system size, the number of modes at ( 7 r, 0) and their relative 
weights are strongly modified by increasing the number of sites. This behavior is consistent 
with the development of a continuum of fractional quasiparticles at ( 7 r, 0) in the thermodynamic 
limit. 

Having established that our Gutzwiller approach depending on wave-vector produces re¬ 
spectively sharp and continuum-like excitations from the |SF) state, we analyze their real-space 
structure to gain further insight into their nature. We consider their overlap with projected 
real-space particle-hole excitations |q, r, +) = P G e* q ' R c^ +rt c R ^|?/ MF ), where a Fourier 
transformation was applied to reflect translation invariance. In this formalism, the most local 
projected particle-hole pair is the spin-flip state S'+|SF) = |q, 0, +) corresponding to a magnon 
while non-local pairs are characterized by a finite separation r. Therefore, the degree of decon¬ 
finement of a fractional S'= 1/2 quasiparticles pair can be characterized using the spatial extent 
of its overlap with projected real-space particle-hole excitations. To evaluate it, we define the 
weighted average 

Pq( r ) = |( f T r ’+l f T ? T+) <q,n,+|S+|S F >| 2 , (6) 

n 

where the aforementioned overlap, represented by the first term, is weighted by the inten¬ 
sity of the transverse spin excitation in the dynamic structure factor. The spatial profile of 
p q (r) for the magnetic zone-boundary wave-vectors, shown in Figs. [l]b and[Tjl, reveals much 
more extended fractional S = 1/2 quasiparticles pairs at (7r,0) than at (7 t/2, 7t/ 2). This is 
confirmed by the system-size dependence of the radially-integrated normalized distribution 
p q( r ) = E|r'|<rPq( r, )> Pitted in Figs. & and[5jl. At (7 t/2, 7t/ 2), P q (r) saturates at a dis- 
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Figure 5: Finite-size effects in the |SF) state. Transverse dynamic spin structure factor at ( 7 r, 0) 
(a) and (7 t/2, 7t/2) (b) for different system sizes ranging from 8x8 (dark blue line) to 32 x 32 
(dark red line). Disk-integrated fractional-quasiparticle pair separation distribution P q (r) at 
(7r, 0) (c) and (7 t/2, 7t/2) (d) for corresponding system sizes, (e) Mean fractional-quasiparticle 
pair separation r q at (zr, 0) (red symbols) and (zr/2, tt/ 2) (blue symbols). The error bars result 
from the variational Monte Carlo sampling. 


tance, r, that is nearly independent of the system size, while at (7 r, 0) it does so at a distance that 
increases with the number or sites. Similarly, the “root-mean-square” fractional quasiparticles 
pair separation r q = |)T/. |r| 2 p q (r)/^ r p q (r)]^ 2 , presented in Fig. JiJe, grows linearly with the 
system size for (7r, 0) while it has a much weaker size dependence for (7t/2, 7t/2). 

Taken together, our real-space results for the | SF) state show that spin excitations at (7r/2, 7 t/ 2) 
can indeed be considered as bound pairs of S' = 1/2 quasiparticles with confined spatial extent. 
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In contrast at ( 7 r, 0) the strong system size dependance of the spin excitations spatial extent 
indicates quasiparticle deconfinement in two spatial dimensions. We do not attempt to extract 
power laws from the numerical data, since the variational |SF) state mimics the long-distance 
spin correlations only qualitatively. Similarly, due to the difficulty in recovering the correct 
asymptotic behavior for the equal-time correlations in the |SF+N) state, it remains an open the¬ 
oretical challenge to obtain explicit quasiparticle deconfinement out of the magnetically ordered 
ground state of Eq. |T[ 

In conclusion, we have presented experimental and theoretical arguments in favor of the de¬ 
confinement of S — 1 excitations into fractional 5 = 1/2 quasiparticles in the excitation spectrum 
of the square-lattice quantum Heisenberg antiferromagnet. Within our Neel plus Staggered-Flux 
ansatz, |SF+N), we show that magnons can be interpreted as bound-pairs of fractional quasipar¬ 
ticles. To the best of our knowledge, the results of our Gutzwiller projection approach provide 
the closest description to date for the high-energy spin-dynamics of the model, as observed 
experimentally in CFTD and relevant for La 2 Cu0 4 . Our polarized neutron scattering results 
unraveled a spin-isotropic excitation continuum at the ( 77 , 0) wave-vector of the square-lattice 
Brillouin zone. Similarily, variational Monte-Carlo calculations in the Staggered-Flux quantum 
spin-liquid ansatz, |SF), reveal an isotropic continuum arising from the spatial deconfinement 
of fractional quasiparticles pairs. While our theoretical findings alone do not unequivocally 
prove the presence of such deconfinement for the exact ground-state of the model, we argue 
that combined with our experimental results, they provide strong support for the presence of 
the analogue of deconfined ID spinons in a two-dimensional non-frustrated long-range ordered 
antiferromagnet. 


15 


References 


[1] Laughlin, R. B. Nobel lecture: Fractional quantization. Rev. Mod. Phys. 71, 863-874 
(1999). 

[2] Su, W. R, Schrieffer, J. R. & Heeger, A. J. Solitons in polyacetylene. Physical Review 
Letters 42, 1698-1701 (1979). 

[3] Hou, C.-Y., Chamon, C. & Mudry, C. Electron fractionalization in two-dimensional 
graphenelike structures. Physical Review Letters 98, 1698-1701 (1979). 

[4] Simon, J. et al. Quantum simulation of antiferromagnetic spin chains in an optical lattice. 
Nature 472, 307-312 (2011). 

[5] Baskaran, G., Zou, Z. & Anderson, P. W. The resonating valence bond state and high-tc 
superconductivity a mean field theory. Solid State Communications 63, 973-976 (1987). 

[6] Balents, L. Spin liquids in frustrated magnets. Nature 464, 199-208 (2010). 

[7] Bethe, H. Zur theorie der metalle. Zeitschrift fiir Physik A Hadrons and Nuclei 71, 205- 
226 (1931). 

[8] Faddeev, L. & Takhtajan, L. What is the spin of a spin wave? Physics Letters A 85, 
375-377 (1981). 

[9] Muller, G., Thomas, H., Beck, H. & Bonner, J. C. Quantum spin dynamics of the an¬ 
tiferromagnetic linear chain in zero and nonzero magnetic field. Physical Review B 24, 
1429-1467 (1981). 


16 



[10] Tennant, D. A., Cowley, R. A., Nagler, S. E. & Tsvelik, A. M. Measurement of the 
spin-excitation continuum in one-dimensional KCuF_3 using neutron scattering. Physical 
Review B 52, 13368-13380 (1995). 

[11] Lake, B., Tennant, D. A., Frost, C. D. & Nagler, S. E. Quantum criticality and universal 
scaling of a quantum antiferromagnet. Nature Materials 4, 329-334 (2005). 

[12] Mourigal, M. el al. Fractional spinon excitations in the quantum heisenberg antiferromag¬ 
netic chain. Nature Physics 9, 435-441 (2013). 

[13] Coidea, R., Tennant, D. A., Tsvelik, A. M. & Tylczynski, Z. Experimental realization of 
a 2D fractional quantum spin liquid. Physical Review Letters 86, 1335-1338 (2001). 

[14] Han, T.-H. et al. Fractionalized excitations in the spin-liquid state of a kagome-lattice 
antiferromagnet. Nature 492 (2012). 

[15] Jeong, M. et al. Field-induced freezing of a quantum spin liquid on the kagome lattice. 
Physical Review Letters 107, 237201 (2011). 

[16] Kozlenko, D. R et al. From quantum disorder to magnetic order in an s=l/2 kagome 
lattice: A structural and magnetic study of herbertsmithite at high pressure. Physical 
Review Letters 108, 187207 (2012). 

[17] Manousakis, E. The spin-| heisenberg antiferromagnet on a square lattice and its applica¬ 
tion to the cuprous oxides. Reviews of Modern Physics 63, 1-62 (1991). 

[18] Anderson, P. W. An approximate quantum theory of the antiferromagnetic ground state. 
Physical Review 86, 694-701 (1952). 

[19] Kubo, R. The spin-wave theory of antiferromagnetics. Physical Review 87, 568-580 
(1952). 


17 



[20] Dyson, F. J. General theory of spin-wave interactions. Physical Review 102, 1217-1230 
(1956). 

[21] Oguchi, T. Theory of spin-wave interactions in ferro- and antiferromagnetism. Phys. Rev. 
117, 117-123 (1960). 

[22] Harris, A. B., Kumar, D., Halperin, B. I. & Hohenberg, P. C. Dynamics of an antiferro- 
magnet at low temperatures: Spin-wave damping and hydrodynamics. Physical Review B 
3, 961-1024(1971). 

[23] Canali, C. M. & Wallin, M. Spin-spin correlation functions for the square-lattice heisen- 
berg antiferromagnet at zero temperature. Phys. Rev. B 48, 3264-3280 (1993). 

[24] Reger, J. D. & Young, A. P. Monte carlo simulations of the spin-1/2 heisenberg antiferro¬ 
magnet on a square lattice. Physical Review B 37, 5978-5981 (1988). 

[25] Hamer, C. J., Weihong, Z. & Arndt, P. Third-order spin-wave theory for the heisenberg 
antiferromagnet. Physical Review B 46, 6276-6292 (1992). 

[26] Anderson, P. W., Baskaran, G., Zou, Z. & Hsu, T. Resonating “ valence-bond theory of 
phase transitions and superconductivity in la 2 cuo 4 -based compounds. Phys. Rev. Lett. 58, 
2790-2793 (1987). 

[27] Auerbach, A. & Arovas, D. P. Spin dynamics in the square-lattice antiferromagnet. Phys¬ 
ical Review Letters 61, 617-620 (1988). 

[28] Hsu, T. C. Spin waves in the flux-phase description of the S=l/2 heisenberg antiferromag¬ 
net. Phys. Rev. B 41, 11379-11387 (1990). 


18 



[29] Ho, C.-M., Muthukumar, V. N., Ogata, M. & Anderson, P. W. Nature of spin excitations in 
two-dimensional mott insulators: undoped cuprates and other materials. Phys. Rev. Lett. 
86, 1626-1629 (2001). 

[30] Coidea, R. et al. Spin waves and electronic interactions in la 2 cuo 4 . Physical Review 
Letters 86, 5377 (2001). 

[31] Headings, N. S., Hayden, S. M., Coidea, R. & Perring, T. G. Anomalous high-energy spin 
excitations in the high-t c superconductor-parent antiferromagnet la 2 cuo 4 . Physical Review 
Letters 105, 247001 (2010). 

[32] Rpnnow, H. M. et al. Spin dynamics of the 2d spin | quantum antiferromagnet copper 
deuteroformate tetradeuterate (cftd). Phys. Rev. Lett. 87, 037202 (2001). 

[33] Christensen, N. B. et al. Quantum dynamics and entanglement of spins on a square lattice. 
Proceedings of the National Academy of Sciences 104, 15264-15269 (2007). 

[34] Singh, R. R. P. & Gelfand, M. P. Spin-wave excitation spectra and spectral weights in 
square lattice antiferromagnets. Physical Review B 52, R15695-R15698 (1995). 

[35] Zheng, W., Oitmaa, J. & Hamer, C. J. Series studies of the spin-! heisenberg antifer¬ 
romagnet at t — 0: Magnon dispersion and structure factors. Phys. Rev. B 71, 184440 
(2005). 

[36] Syljuasen, O. F. & Rpnnow, H. M. Quantum renormalization of high-energy excitations 
in the 2D heisenberg model. Journal of Physics: Condensed Matter 12, L405 (2000). 

[37] Sandvik, A. W. & Singh, R. R. P. High-energy magnon dispersion and multimagnon 
continuum in the two-dimensional heisenberg antiferromagnet. Phys. Rev. Lett. 86, 528- 
531 (2001). 


19 



[38] Liischer, A. & Lauchli, A. M. Exact diagonalization study of the antiferromagnetic spin- 
1/2 heisenberg model on the square lattice in a magnetic field. Physical Review B 79, 
195102 (2009). 

[39] Zhitomirsky, M. E. & Chernyshev, A. L. Colloquium : Spontaneous magnon decays. Rev. 
Mod. Phys. 85, 219-242 (2013). 

[40] Tang, Y. & Sandvik, A. W. Confinement and deconfinement of spinons in two dimensions. 
Phys. Rev. Lett. 110, 217213 (2013). 

[41] Kohno, M., Starykh, O. A. & Balents, L. Spinons and triplons in spatially anisotropic 
frustrated antiferromagnets. Nature Physics 3, 790-795 (2007). 

[42] Enderle, M. et al. Two-spinon and four-spinon continuum in a frustrated ferromagnetic 
spin-1/2 chain. Physical Review Letters 104, 237207 (2010). 

[43] Gros, C. Superconductivity in correlated wave functions. Phys. Rev. B 38, 931-934 (1988). 

[44] Gros, C. Physics of projected wavefunctions. Annals of Physics 189, 53 - 88 (1989). 

[45] Dmitriev, D. V., Krivnov, V. Y., Likhachev, V. N. & Ovchinnikov, A. A. Variation function 
with vortexes in the heisenberg 2-dimensional anti-ferromagnetic model. Phys. Solid State 
38, 397 (1996). [Translated from: Fiz. Tv. Tela 38, 397 (1996)]. 

[46] Wen, X.-G. & Lee, P. A. Theory of underdoped cuprates. Phys. Rev. Lett. 76, 503-506 
(1996). 

[47] Nayak, C. Density-wave states of nonzero angular momentum. Phys. Rev. B 62, 4880- 
4889 (2000). 


20 



[48] Lee, T. K. & Feng, S. Doping dependence of antiferromagnetism in la 2 cuo 4 : A numerical 
study based on a resonating-valence-bond state. Phys. Rev. B 38, 11809-11812 (1988). 

[49] Trivedi, N. & Ceperley, D. M. Ground-state correlations of quantum antiferromagnets: A 
green-function monte carlo study. Phys. Rev. B 41, 4552-4569 (1990). 

[50] Calandra Buonaura, M. & Sorella, S. Numerical study of the two-dimensional heisenberg 
model using a green function monte carlo technique with a fixed number of walkers. Phys. 
Rev. B 57, 11446-11456 (1998). 

[51] Paramekanti, A., Randeria, M. & Trivedi, N. High-T c superconductors: a variational 
theory of the superconducting state. Phys. Rev. B 70, 054504 (2004). 

[52] Ivanov, D. A. Resonating-valence-bond structure of gutzwiller-projected superconducting 
wave functions. Phys. Rev. B 74, 024525 (2006). 

[53] Syromyatnikov, A. V. Spectrum of short-wavelength magnons in a two-dimensional quan¬ 
tum heisenberg antiferromagnet on a square lattice: third-order expansion in 1/s. Journal 
of Physics: Condensed Matter 22, 216003 (2010). 

[54] Li, T. & Yang, F. Variational study of the neutron resonance mode in the cuprate super¬ 
conductors. Phys. Rev. B 81, 214509 (2010). 

Acknowledgments 

We gratefully acknowledge fruitful discussions with M. Zhitomirsky, S. Sachdev, C. Broholm 

and L. P. Regnault. Work in EPFL was supported by the Swiss National Science Foundation, 

the MPBH network, and European Research Council grant CONQUEST. Computational work 

was supported by the Swiss National Supercomputing Center (CSCS) under project ID s347. 


21 



Work at Johns Hopkins University was supported in part by the US Department of Energy, 
office of Basic Energy Sciences, Division of Material Sciences and Engineering under grant 
DE-FG02-08ER46544. NBC was supported by the Danish Agency for Science, Technology 
and Innovation under DANSCATT. 


22 



Supplementary information for: Fractional 
excitations in square-lattice quantum 
antiferromagnet 


B. Dalla Piazza, 1 ’* M. Mourigal, 1,2,3 ’* N. B. Christensen, 4,5 
G. J. Nilsen, 1,6 P. Tregenna-Piggott, 5 T. G. Perring, 7 
M. Enderle, 2 D. F. McMorrow, 8 D. A. Ivanov, 9,10 and H. M. R0nnow' 

laboratory for Quantum Magnetism, 

Ecole Polytechnique Federale de Lausanne (EPFL), CH-1015, Switzerland 
2 Institut Laue-Langevin, BP 156, 

F-38042 Grenoble Cedex 9, France 

^Institute for Quantum Matter and Department of Physics and Astronomy, 

Johns Hopkins University, Baltimore, MD 21218, USA 
4 Department of Physics, Technical University of Denmark (DTU), DK-2800 Kgs. Lyngby, Denmark 
’Laboratory for Neutron Scattering, Paul Scherrer Institute, 

CH-5232 Villigen PSI, Switzerland 
department of Chemistry, University of Edinburgh, 

Edinburgh, EH9 3JJ, United Kingdom 
' ISIS Facility, Rutherford Appleton Laboratory, 

Chilton, Didcot, Oxon 0X11 OQX, United Kingdom, 

8 London Centre for Nanotechnology and Department of Physics and Astronomy, 

University College London, London WC1H OAH, United Kingdom 
institute for Theoretical Physics, ETH Zurich, 

CH-8093 Zurich, Switzerland 
10 Institute for Theoretical Physics, University of Zurich, 

CH-8057 Zurich, Switzerland 

*To whom correspondence should be addressed; E-mail: bastien.dallapiazza@epfl.ch, mourigal@jhu.edu 


23 



Contents 


1 Experimental Materials and Methods |25] 

1.1 Properties of the metal-organic compound CFTD. |25] 

1.1.1 Crystal structure . [25] 

1.1.2 Crystallographic notations. |25] 

1.1.3 Magnetic properties . |25] 

1.1.4 Magnetic order. |26] 

1.1.5 Magnetic excitation spectrum . |27] 

1.2 Crystal Growth. |27] 

1.3 Time-of-flight inelastic neutron scattering experiment. |28] 

1.4 Polarized triple-axis inelastic neutron scattering. |29] 

1.4.1 Experimental set-up . |29] 

1.4.2 Polarization analysis. |29] 

1.4.3 Reconstruction of the dynamic structure factor. |30] 

1.4.4 Data correction. |30] 

1.4.5 Results. |32] 

2 Theoretical Materials and Methods |33] 

2.1 SF+N mean-field Hamiltonian. |33] 

2.2 Numerical projection of the Heisenberg Hamiltonian. |34] 

2.3 Transverse dynamic structure factor. |35] 

2.4 Transverse magnon mode intensity. |36] 

2.5 Longitudinal dynamic structure factor. |37] 

2.6 Neel-held effect on fractional quasiparticle deconfinement. |38] 


24 






















1 Experimental Materials and Methods 

1.1 Properties of the metal-organic compound CFTD 

A detailed review of the properties of CFTD can be found in Chapter 3 of Ref. 0521. 

1.1.1 Crystal structure 

The layered metal-organic compound Cu(DC00) 2 -4D 2 0 (CFTD) investigated in this work is 
the deuterated analog (hydrogen 1 H substituted by deuterium 2 D) of copper formate tetrahydrate 
(CFTH) 0531 . It is described in the monoclinic space-group P2!/a at 300 K and in P2 i/n below 
T s « 246 K where the unit-cell doubles along c 0541 . This change of symmetry results from a 
first-order antiferroelectric ordering involving the water molecules between the copper-formate 
planes, see Fig. While deuteration influences the temperature of the structural transition 
T s , CFTH and CFTD are structurally very similar 0541 . At T = 120 K, the lattice parameters 
of CFTD are a = 8.113 A, b = 8.119 A, c = 12.45 A and (3 = 100°. The two non-equivalent 
copper sites in the unit-cell are coordinated to four oxygens from formate groups in the ab plane 
and to two oxygens from crystalline water above and below that plane, forming a staggered 
pattern of elongated octahedrals. The elongated directions (labeled ||) are within « 25° of 
c. The two-dimensional lattice of Cu 2+ ions is face centered with nearest-neighbor coppers 
arranged on a square lattice (within 0.07%) at an average distance d = \\Ja? + b 2 ~ 5.74 A, 
see Fig. €)>• The monoclinic angle j3 influences the stacking of the ab planes along c but is 
irrelevant for the two-dimensional physics considered in this work. Assumming the two copper 
sites indistinguishable, we define the elementary two-dimensional square-lattice by the nearest 
neighbor vectors x = (a + b)/2d and y = (a — b)/2d such that |a;| = \ y\ = 1. 

1.1.2 Crystallographic notations 

In the reciprocal lattice of CFTD spanned by a*, b*, and c*, we define wave-vectors as Q = 
ha* + kb* + ic* with h, k and l in reciprocal lattice units (r.l.u). In the i = 0 plane considered 
in this work, Brillouin zone centers (Q = r) are defined by the equivalent reciprocal-lattice 
vectors r = (2,0,0), r = (1,1,0) and r = (0,—2,0). In turn, we define the elementary 
reciprocal square-lattice spanned by x* and y* such that the two-dimensional wave-vector q = 
q x x* + q y y* with q x = n(h + k) and q y = ir(h — k) conform with standard theoretical notations. 

1.1.3 Magnetic properties 

The magnetic susceptibility of CFTH/CFTD reveals an overall antiferromagnetic behavior with 
a negative Weiss constant 0 W = — 175 K 055ll . a 2D short-range order indicated by a broad maxi- 
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Supplementary Figure 1: Crystal and magnetic structure of CFTD. (a) Three dimensional 
crystal structure with Cu, O, C and H/D atoms in blue, red, brown, gray, respectively. Note that 
the formate groups in the z = 0 and z — 1 planes are not represented, (b) Two-dimensional 
copper-formate plane in CFTD and dominant superexchange path J between copper ions, (c) 
Sketch of the magnetic structure of CFTD. Note that the monoclinic angle and the out-of-plane 
spin direction have been amplified for clarity. 

mum around T* « 65 K Ii57ll . and a transition to 3D long-range order at T N = 16.5(5) K If56ll . The 
analysis of the T > 35 K susceptibility using high-temperature series-expansion yields a nearest- 
neighbor exchange J = 6.2(3) meV ll58l . Results from local probes such as NMR Ii59l, 60], 
ESR Ml and isothermal magnetization ll62l [63 1 conform with this picture. To fully account 
for these measurements, it is however necessary to introduce small gyromagnetic and exchange 
anisotropies. The gyromagnetic-tensor of the Cu 2+ sites is staggered, with estimated princi¬ 
pal components g± = 2.1 and g\\ = 2.4, and average value of g av = 2.19 [ 571 [621. In turn, 
the dominant nearest-neighbor exchange J is weakly anisotropic, with an estimated symmetric 
off-diagonal component < 10~ 3 J ll611 and an estimated antisymmetric component < 1(1 2 J 
with Dzyaloshinskii-Moriya vector directed in the ac-plane, closer to c than to a [f62l 163 1. The 
interlayer exchange is estimated to be as small as J c « 5 x 10 ~ 5 J. Overall, CFTD is one of 
the best known realization of the spin- 1/2 square-lattice Heisenberg antiferromagnet, with an 
exchange energy-scale perfectly suited for polarized neutron scattering experiment, in contrast 
to the much higher energy-scale of the parent compounds of the superconducting cuprates. 

1.1.4 Magnetic order 

Unpolarized and polarized neutron diffraction measurements 11641 establish that CFTD indeed 
hosts long-range magnetic ordering below T N . The static magnetic structure is best described 
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by an anticollinear arrangement of nearest neighbor spins, with spins lying in the ac-plane with 
an on-site moment (m) of norm 0.48(2) and direction 8(1)° away from a towards c, i.e. 
only 3(1)° from a*. The four-sublattice structure, see Fig. allows a very weak canting 
along b but a net ferromagnetic contribution > 0.005(6) fi B is excluded. This suggests that 
the magnetic structure is stabilized by a small off-diagonal symmetric exchange with negligible 
antisymmetric exchange. Only magnetic reflections with £ = 2n + 1 are observed so that there 
is no magnetic Bragg peak accessible in the a*b* reciprocal plane. The lowest-angle magnetic 
reflection is Q = (1,0,1). As the ^-component of Q is irrelevant for the two-dimensional 
correlations considered in this work, Q = (1, 0, 0) and equivalent reflections correspond to the 
center of magnetic Brillouin zone (M-point) with q = (7r, 7 r). Similarly, Q = (0.5, 0.5, 0) and 
equivalent positions correspond to the X-point of the magnetic zone-boundary point q = (M), 
while Q = (0.5, 0, 0) and equivalent positions coincide with the zone-boundary point q = 
(tt/2, tt/2). 

1.1.5 Magnetic excitation spectrum 

Previous inelastic neutron scattering studies on single-crystals of CFTD confirmed it is an ex¬ 
cellent realization of the Heisenberg square-lattice antiferromagnet [|65] [66, 67; 68, 69 ]. A fit to 
the dispersion of its magnetic excitations using spin-wave theory, including the renormalization 
factor Z c = 1.18, yields J = 6.11(2) meV [|69|| , in good agreement with susceptibility results. 
Likewise, exchange terms going beyond that of the spin-1/2 square-lattice Heisenberg antifer¬ 
romagnet are found to be very small. No dispersion is observed along c* 11651 what confirms 
the excellent two-dimensionality of the coumpound. The largest deviation from the ideal model 
is a small gap at the magnetic-zone center Q = (1, 0, 0) of « 0.38 meV [66,67', [68], attributed 
to a small off-diagonal exchange interaction. The absence of further-neighbor interactions is 
inferred indirectly, by comparing the observed dispersion at the magnetic zone-boundary with 
theoretical results Ii67ll69l 

1.2 Crystal Growth 

Due to the large incoherent neutron cross section of 1 H. it was necessary to prepare a 2 D- 
substituted sample of CFTD for our neutron scattering experiments. This was achieved by 
dissolving CuO in a solution of D 2 -formic acid in D 2 0 and evaporating the solution to pre¬ 
cipitate CFTD. The product was purified by redissolving, filtering, and recrystallizing twice. 
Seed crystals were grown by slow evaporation of a supersaturated solution from beakers coated 
with dime thy ldichloro silane. All operations were carried out under inert atmosphere to avoid 
exchange with airborne water. 

The large single crystals used in the experiments were grown by a convection method ifTPl : 
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two glass columns were charged with deuterated starting material and a seed crystal suspended 
on a glass fiber, respectively, and connected at top and bottom to form a loop. The apparatus 
was filled with a saturated solution of CFTD in D 2 0, and a convection current was induced 
by regulating the temperatures of the starting material column with a water bath. The flow of 
supersaturated solution generated in the warmer starting material zone over the seed crystal in 
the cooler zone lead to the growth of a 7.2 g cubic crystal of dimension 20 x 20 x 20 mm 3 
within a period of around 14 days. Three smaller ~ 4 g crystals were obtained in subsequent 
growth cycles. 

1.3 Time-of-flight inelastic neutron scattering experiment 

Our time-of-flight inelastic neutron scattering results were obtained using the chopper spec¬ 
trometer MAPS at the ISIS Facility, Rutherford Appleton Laboratory (UK) with an incident 
neutron energy E, = 36.8 meV. The sample consisted of three co-aligned crystals of total mass 
~ 12 g mounted on an aluminum sample holder, aligned with the reciprocal directions a* and 
b* kept perpendicular to the direction of the incident beam, and cooled down to T = 6 K in 
a closed-cycle cryostat. After standard corrections and transformations, the intensity obtained 
at momentum transfer /iQ and energy transfer two is proportional to diagonal components of 
the dynamic structure factor S aa ( Q,co) defined as the time Fourier-transform of the thermally 
averaged spin-spin correlation function, 

i r°° 

S aa (Q,co) = — J dte~ tlJJt (S a (Q, 0)S a (Q,t)) T , (7) 

with S^Q, t) = l/(27r) 3 f S a (r, t) exp(?'Q • r)d 3 r the space Fourier transform of the a com¬ 
ponent of the spin operator S(r,t). In our time-of-flight experiment a« components that are 
perpendicular to Q contribute to the intensity. As an array of detectors is used to cover a large 
solid angle, each pixel records a different linear combination of oa-compomcnts. 

Benefiting from the absence of dispersion of the excitations along c*, the collected data 
were projected in the a* b* reciprocal-plane and integrated along the third momentum transfer 
direction. Symmetries of the elementary square-lattice were used to fold the data into an ir¬ 
reducible 1/8 of the two-dimensional Brillouin zone. An energy-dependent background was 
subtracted from the obtained data. Originating from the imperfect deuteration of our sam¬ 
ple due to 1 H/ 2 D exchange, this coherent and nuclear spin incoherent phonon background 
is Q independent within a good approximation and obtained by integrating on a square Ah, 
A k = ±0.125 r.l.u around the nuclear zone center r = (1,1, 0). 
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1.4 Polarized triple-axis inelastic neutron scattering 

1.4.1 Experimental set-up 

Our polarized inelastic neutron scattering results were obtained using the thermal-triple-axis 
neutron spectrometer IN20 at the Institut Laue-Langevin, Grenoble (France) equipped to per¬ 
form longitudinal xyz polarization analysis. The sample was mounted on an aluminum sample 
holder with the reciprocal axes a* and b* in the scattering plane and inserted in a standard 4 He 
cryostat reaching base temperature T = 1.5 K. Polarized neutrons were produced and analysed 
by horizontally focused Heusler (111) monochromator and analyzer crystals with fixed vertical 
curvature. The final neutron wave-vector was fixed at kf = 2.662 A -1 . No collimators were 
installed and the instrument was operated in W configuration. A spin flipper was installed in 
the scattered beam along with a PG filter to suppress second order contamination. A neutron 
flux-monitor was installed in Ay to normalize the collected intensity. 

1.4.2 Polarization analysis 

In the above Heusler-Heusler configuration, a set of Helmholtz coils around the sample allowed 
to provide a guide field of 10-60 Gauss at the sample position, and thus to define the direction 
of the neutron polarization at the sample. The field and polarization direction is chosen either 
parallel to the momentum transfer Q (x {) ), perpendicular to Q within the scattering plane (y 0 ) or 
perpendicular to the scattering plane, (z (l ). Scattering without a change of the neutron spin state 
(non-spin flip scattering) is denoted accordingly with labels x 0 x 0 , yoyo, or z 0 z 0 , respectively, 
while spin-flip scattering is measured by activating a neutron spin-flipper in the outgoing beam 
after the sample, and labeled x 0 x^, yoyo, z 0 Zo- 

This way, we measured the following partial cross sections: 
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where {A* A) 

7> 

denotes the time Fourier transformed and thermally averaged pair correlation 


function, 

(A*A) t „ = — / dte~**(A*(Q,0)A(Q,t)) T , (9) 

J —oo 
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N( Q) is the nuclear scattering amplitude, -Mj_q(Q) the magnetic scattering amplitude per¬ 
pendicular to the momentum transfer Q, u iso denotes the nuclear spin-independent elastic and 
inelastic incoherent scattering due to random isotope distribution, and <j nsi the nuclear spin- 
dependent elastic and inelastic incoherent scattering related to random nuclear spin orientation. 
The proportionality constant is the same for all six equations. From the above six equations, 
the four quantities (cr ns i}T,co, and (N*N) T}U} +(u iso ) T:U can be 

determined separately. 

1.4.3 Reconstruction of the dynamic structure factor 

The components of the magnetic vector M° (Q. t) entering the polarization analysis formulas of 
Eq. Sj8]are related to the space Fourier-transform of spin operators, M ,y (Q. t) = ^gf(Q)S a (Q, t ), 
where g is the the gyromagnetic factor and /(Q) the magnetic form-factor of the considered 
Cu 2+ ions. Polarization analysis thus allows to extract separately two components of the dy¬ 
namic structure factor at each Q. 

Note that the a = x,y,z components of 5“"(Q,o;) are expressed in the Cartesian co¬ 
ordinate system attached to the ordered-moment direction 2 || (m) while the 6 = yo . zq 
components of M_q(Q) are expressed in the Cartesian coordinate system attached to the 
momentum transfer .x 0 || Q. This allows to reconstruct the full dynamic structure factor 
by combining measurements from two equivalent Brillouin Zones, the orientations of which 
with respect to (m) differs, see Fig. S[2j In the r = (2, 0, 0) zone, Q is mostly parallel to 
(m) what allows to obtain the transverse components S xx (Q. u) and S yv (Q. u) by measuring 
Qi = (2.5, 0, 0) and Q 2 = (2.5, —0.5, 0) corresponding to the 2D wave-vectors q = (n/2, n/2) 
and q = (7r, 0), respectively. Fikewise, in the r = (0, —2,0) zone, Q is mostly perpendic¬ 
ular to (m), what allows to obtain the longitudinal and transverse components S ZZ (Q, oS) and 
<S TO (Q, oj) for Q 3 = (0, —2.5, 0) and Q 4 = (0.5, —2.5, 0) corresponding to the 2D wave-vectors 
q = (7 t/2,7t/ 2) and q = (7r,0), respectively. This is conveniently understood schematically, 
see Fig. S(2j 

1.4.4 Data correction 

Our experiment thus consisted in collecting the six above partial cross sections of Eq. Sj8] as 
a function of energy transfer for the four momentum transfers of Fig. S0 Various corrections 
had to be applied to relate the measured intensities to the components of the dynamic structure 
factor. First, the intrinsic k^ 1 sensitivity of the flux monitor, corrected for higher-order neutron 
contamination, balanced the /c* dependence of the neutron scattering cross-section. In addition, 
a geometrical correction was applied to account for the variation in the vertical focal length at 
the sample position as a function of incident neutron energy. Then, due to the imperfect beam 
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Supplementary Figure 2: Representation of the a*b* reciprocal plane of CFTD. The gray areas 
represent the nuclear zones (center T) while the white areas are the magnetic zones (center M). 
The magnetic zone-boundaries are represented as thick solid (red and blue) squares. The wave- 
vectors Q investigated in our triple-axis experiment are represented by green circles along with 
their corresponding 2D wave-vector q from the closest nuclear zone center r. The direction of 
the ordered moment (3(1)° from a* towards c) is represented as a bold arrow. The “transverse” 
and “longitudinal” configurations with respect to the ordered moment are schematically drawn 
in blue and red, respectively. 


polarization, the data were corrected for the mutual influence of non-spin-flip and spin-flip 
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scattering into another in each channel, 
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where Pg is the beam polarization in the direction (3, d 2 a/dfldEf is a measured cross-section 
and d 2 a/dPldEf\' a corrected one. Following, careful tests performed during the experiment 
on the polarization of the direct-beam, we used an isotropic imperfect polarization P XQ = P yo = 
P Z0 =0.873(5). Given the small anisotropy of the gyromagnetic factor and noting that all wave- 
vectors considered here have similar length within 2%, we obtain that the magnetic form factor 
g f (QO/2 is constant. We also neglected the difference in extinction between our four Q, given 
the nearly cubic shape of the sample and also neglected the 3° out-of-plane orientation of the 
ordered moment. 


1.4.5 Results 

Finally, we obtain (M^qM|° q ) t (Q, u>) and (M| 1 qM|°q) t (Q, cc) after correcting the measured 
cross-sections using Eq. Sj8j and reconstruct the full dynamic structure factor using the equations 
below. For the q = (n/2, 7t/ 2) wave-vector, we reconstruct the dynamic structure factor using: 

5“(q,w) = 

S”(q,w) = i[<^M5,) T (Q 1 ,u) + <M^MS,) r (Q 3 ,u)]. 

For the q = (n, 0) wave-vector, we reconstruct the dynamic structure factor using 

S Z2 (q,w) \ _ / sin 2 6 a cos 2 £ a \ 1 / (M^Mf Q ) T (Q 2 , u) \ 

S xx (q,u) ) \ cos 2 6b sin 2 0 b J ^ (MIqM^° q ) t (Q 4 ,uj) J 

S m (q,u,) = j[(M^MS,) t (Q 2 ,w) + (M^M5,) t (Q 4 ,u)] , 

where 9 a m 9 b ^ 10° accounts for the angle between Q and the directions of the reciprocal space 
vectors a* and b*, respectively. In the end, the transverse, longitudinal and total contributions 
to the dynamic structure factor are obtained as, 

^(qjw) = ^[«S +- (q,o;) +<S~+(q,o;)] = <S X3; (q, w) + S ra (q,u;), 

5 L (q,a;) = S zz (q,u), (13) 

5(q,w) = 5**(q,a;)+5*®(q,a;)+5«(q,a;). 


(ID 


( 12 ) 
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2 Theoretical Materials and Methods 


2.1 SF+N mean-field Hamiltonian 


The SF+N Hamiltonian is a mean-field ansatz for the Heisenberg model expressed in terms of 
fermionic operators (Eq. [2]). The mean-field decoupling is: 

■HmF = (Xv'4r C 3° + h - C ) + ^{- l ) l Kc\ a C i!T (14) 

ia 


with 


and 


h a = -<jH n 



Xij = e i6ij dij = (-iy* +j y0 o 


(15) 

(16) 


for nearest-neighbor pairs (i,j), i x ( y ) being the x (y) coordinate of site i on the square lattice. 
For constructing the ground-state wave function, we ignore the overall scale of the Hamiltonian 
"Hmf and put \xij\ = 1. The two variational parameters are 8 0 and fF N (the Neel field). The 
complex hopping amplitude Xij causes the square lattice to be threaded with staggered fluxes 
going through the lattice squares. Electrons circulating around a square acquire a phase of ±400. 

The mean-field Hamiltonian "Hmf is diagonalized by the quasiparticle operators 


Tkir— 

Tkcr+ 


^kcr— 


u 


kcr+ 



(^kcr “h Ck+na) ^k a— 
(^kcr “1“ Ck+ncr) “t“ ^kcr+ 



(Qccr 
(Qc a 


Ck+ncr) j 
Ck+ncr) ? 


(17) 


where the wave vector k is restricted to the magnetic Brillouin zone, a is the spin index, and 
n = ( 7 r, 7r). The coefficients u krj± and v k(T ± are calculated as 


^kcr— ^k <7+ 


^kcr— — ^kcr+ — 



(18) 


Here ± subscripts correspond to the upper and lower branch of quasiparticles with energies 
±cu k , respectively. Note that these energies do not carry any physical meaning, but only deter¬ 
mine the empty and occupied quasiparticle states. They are given by 

cu k = \J |A k | 2 ± , A k = 1 ( e l6 ° cos k x ± e~ l6 ° cos k y ) . (19) 
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The variational parameters 0 o and //\ are optimized using the standard variational Monte Carlo 
technique. The optimal wavefunction |SF + N) is obtained with 

6*0 = 0.17T 

Hx = 0.11 

On the other hand, the spin liquid state |SF) is obtained by letting the Neel field go to zero 

6 0 = O.Itt (22) 

H n = 0. (23) 


( 20 ) 

( 21 ) 


2.2 Numerical projection of the Heisenberg Hamiltonian 


In our variational approach, we project the physical Heisenberg spin Hamiltonian T~L (given by 
Eq. [Tj) onto the subspace spanned by the Gutzwiller-projected mean-field particle-hole states 
Ik, q) (Eq. For simplicity, in this section we only discuss the case of the transverse excita¬ 

tions; the longitudinal case is treated in a similar way. 

Diagonalizing the projected Hamiltonian amounts to solving a generalized eigenvalue prob¬ 
lem 

E «&<'«', = E m E < 24 > 

k' k' 

in order to determine the coefficients ©F of the expansion of the excitations in the basis | k, q) 
(Eq. |4j) and the excitation energies E qn . The matrices T~L kk , and O kk , are defined as 

^kk' = ( k > q|^l k/ ; q> > c’kk' = ( k > ql k '> q) • (25) 


These matrices are sampled using the Monte Carlo reweighing technique developed by Li and 
Yangf54j|: 


n 


q 

kk' 


o 


q 

kk' 


Tf( nq) v K k ’ q l a )l 2 v ( k -ql a )( Q 'l^l/ 3 )(/ j l k/ ^q) 
1 J ^jrEk" l( k ">q|«)l 2 

Pq(“) 


Ti(O q )E^( Q ) 

a 


(k,q|a)(a|k',q) 

Ek"l( k "^q|«)| 2 ‘ 


(26) 


Here |a) and | (3) are real space spin configurations and p q (a) are the Monte Carlo sampling 
probabilities. The amplitudes (k, q|a) are calculated as Slater determinants [l44ll . The ampli¬ 
tudes (/31k', q) are efficiently calculated by changing a row and a column in the Slater matrix 
using a rank-2 determinant update. Note that this approach leaves uncalculated the overall nor¬ 
malization factor Tr(C> q ) = V k O kk in Eq. S26 which can be ignored for the generalized 
eigenvalue problem (Eq. 
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2.3 Transverse dynamic structure factor 


Since the local spin-flip operator S'/ commutes with the Gutzwiller projection, its action on 
the variational ground state |SF(+N)) can be decomposed in the projected particle-hole states 
|k,q): _ 

SJ|SF(+N)) = £<A kq |k,q), (27) 

k 


with the expansion coefficients 


^kq ^kf+^k—q4_— d~ ^kt+^k—qj,— • 


(28) 


More generally, a projected transverse particle-hole pair state can be written as 

|q,r,+) = P G J^e* q ‘ R cJ l+rt c R 4,|V’ M F) = X^ k n( r )l k > q ) 

R k 


with 


(29) 


0kq(/) 6 /r (^kf+^k— qj,— d - ^kf+^k—qj,—) d - G (^k|+^k—qj,— d - ^kt+^k—qj,—) ] i (30) 

e r = 2 (l + e ?n ‘ r ) , and e r = \{l — e ,n ‘ r ) . The zero-temperature transverse dynamic structure 
factor of a system in the groundstate \ipGs) (with the ground-state energy E GS ) is 

s ± (q,w) = X)(V’Gs|S-|A q ><\|SJIV’GS>«(« - E Xq + E cs ) , (31) 

Aq 


where {| A q )} is the set of all excited states with energies E Xtl - In our variational approxima¬ 
tion, we restrict the sum to the projected particle-hole eigen-states as obtained by solving the 
generalized eigenvalue problem (Eq. S|24]) and obtain: 


^(q.w) 


E 


E' 

k,k' 


'n* /7-)U 
’kq^kk' 1 


Vq 


5(tu - E+ n + E SF ( +N )). 


(32) 


Note that in order to calculate this expression, we need to know the normalization of the overlap 
matrix (T kk ,. This normalization was disregarded in the generalized-eigenvalue calculation (Eq. 


S 26), but can be calculated independently from the sum rule 


{S q S+) = / dw5 ± (q,o;) = 




E' 

k,k' 




k.k' 


■ ( 33 ) 


This is the equal-time transverse spin correlation function, which can be calculated directly as 
an expectation value in the variational ground state by using a Monte Carlo sampling. This 
method allows us to calculate the correctly normalized transverse dynamic structure factor. 
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Supplementary Figure 3: Magnon intensity along the high-symmetry directions of the Brillouin 
zone. Solid black line: linear spin-wave theory. Red dots: Lowest energy exciations of the 
|SF+N) variational state. Blue errorbars: experimental results from ref. Il33ll . The |SF+N) and 
linear spin-wave theory intensities are scaled such that they coincide with the experimental data 

at q = (tt/2, tt/2). 

2.4 Transverse magnon mode intensity 

We compare here our variational calculation on the |SF+N) state with SWT. Remarkably, it 
reproduces well not only the energies of the magnon excitations, but also the intensities 

^q,n = Kq^lSq |V>Gs)| 2 , (34) 

except for the low-energy Goldstone magnons and for the (7r, 0) point (where some part of the 
calculated spectral weight is transferred to higher energy states), see Fig. 50 However both 
theories fail to capture the dramatic loss of spectral weight observed at (7r, 0) in experiments, as 
explained in the main text. 
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Supplementary Figure 4: The longitudinal dynamic structure factor of the |SF+N) state for the 
system size 24 x 24. The dashed white line shows the energies of the magnon mode found in 
the transverse dynamic structure factor. 

2.5 Longitudinal dynamic structure factor 

To obtain the longitudinal dynamic structure factor, we construct a set of A S z = 0 longitudinal 
excited states: 

|k, q, era) = P G 7]L+7k-qa- \^mf) (35) 

At q = (0, 0) and q = ( 7 r, 7r), we must complement this set with the ground state P G |^mf) 
since they all belong to the same momentum in the magnetic Brillouin zone and thus might 
overlap. The numerical projection of the Heisenberg Hamiltonain on this set of states is per¬ 
formed in exactly the same way as in the transverse case. Note that calculating the amplitudes 
(/31 k' , q, a a) now involves updates of the Slater determinant by changing for instance two rows 
and one column simultaneously, which can be done by generalizing the determinant update 
formula to arbitrary rank-n change of rows and columns. 

To calculate S zz (q, ui), we derive the coefficients 0£ qcr such that 

S^mf) = ^ 0kqJ k > a(T ) +^n^ — |^mf) • (36) 

■ ™ ‘ ™ (ih. 
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An explicit calculation gives 

0kq CT = ^ K,+«k-q ff - + ^ka+^’k-qa-) • (37) 

As in the transverse case, the longitudinal dynamic structure factor is obtained by restricting 
the sum over the excited states to the set of projected longitudinal particle-hole states (Eq. S |35| ) 
giving 

S“(q,w) = Yi |(q,n,0|S‘|GS>| 2 ^ - E° v + E os ) (38) 

k 

where |q, n, 0) is the longitudinal eigenstate corresponding the eigenvalue E[ yn obtained by di¬ 
agonalizing the Heisenberg spin Hamiltonian (Eq. |Tj) onto the non-orthonormal states |q, k, a a) 
(Eq. S(35>. 

The normalization of the longitudinal structure factor is performed in the same way as in 
the transverse case from the sum rule 

J M»(q, W ) = (SVq)- (39) 

In Fig. S|4j we show the longitudinal dynamic structure factor calculated in the |SF+N) state. 
Note that for the |SF) state, the longitudinal structure factor coincides with the transverse one. 

2.6 Neel-field effect on fractional quasiparticle deconfinement 

In this section, we provide additional data on the variational excitations and ground-state prop¬ 
erties as a function of the variational field H^. As IIs interpolates between the |SF) state 
(6* 0 = 0.17T, H n = 0) and the |SF+N) state (the optimal values 9 0 = 0.17T and i/ N = 0.11), the 
transverse spin fluctuations change from algebraic to exponential (Fig. Sj5^), the continuum in 
the variational spectrum is pushed to higher energies leaving behind a magnon peak (Fig. Sj5J>), 
and the spatial extent of the S' = 1/2 quasiparticles pairs decreases (Figs. Sg and S(5ji). This 
behavior is consistent with our conjecture that the long-distance transverse spin-fluctuations 
play a crucial role for the possibility of fractional quasiparticle deconfinement. 
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Supplementary Figure 5: Dependence of the ground-state spin correlations and excitations at 
(7r, 0) in the |SF+N) state on the variational Neel field // N -. The system size is 32 x 32 and 
the value of the staggered flux is kept at the optimal level for the |SF+N) state (0 O = 7r/10). 
(a) The equal-time transverse spin correlation function in real space in the log-log scale, (b) 
The transverse dynamic structure factor, (c) The disk-integrated quasiparticle-pair separation 
distribution (defined in the same way as in the panels c and d of Fig. [5]). (d) The root-mean- 
square quasiparticle-pair separation (defined in the same way as in the panel e of Fig. [5]). Lines 
and marker colors corresponds to panel C legend. 
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